1 Especies invasoras - Un problema global
- 1.1 ¿Por qué son invasoras?
2 Objetivos del trabajo y definiciones iniciales
- 2.1 Objetivos
- 2.2 Definiciones
3 Paquetes necesarios
4 Delimitación del MCI
5 Análisis exploratorio
6 Descartar variables muy correlacionadas (opcional)
7 Reducción de dimensionalidad mediante PCA
- 7.1 PCA con todas las variables
- 7.2 PCA con la función ‘VIF’
8 Obtención de datos de países y proyección sobre el MCI
9 Cuantificación del grado de similitud climática
10 Determinación del riesgo invasor
11 Conclusiones
Bibliografía

1 Especies invasoras - Un problema global

El crecimiento poblacional y la industrialización han llevado al hombre a ocupar prácticamente todos los biomas, impactando gravemente en la diversidad de los sistemas naturales. La biodiversidad es esencial para el mantenimiento de las funciones y servicios ecosistémicos, de los cuales nos beneficiamos, pero a la vez son indispensables para el mantenimiento de la homeóstasis terrestre (Cardinale et al. (2012)). Los esfuerzos para la conservación de estas funciones han llevado a analizar las causas de pérdida de la biodiversidad, siendo las especies invasoras la segunda mayor amenaza a la biodiversidad, después de la pérdida de hábitat (Simberloff et al. (2013)). Para luchar contra esta causa necesitamos políticas de gestión efectivas, pero el esfuerzo de erradicación aumenta con el tiempo de invasión, y por este motivo es imprescindible tratar de trabajar en la previsión y primeras fases de invasión (Van Wilgen, Davies, and Richardson (2014)).

Según la Unión Internacional para la Conservación de la naturaleza (UICN), las especies invasoras son:

“Especie exótica invasora”: especie exótica que se establece en un ecosistema o hábitat natural o seminatural; es un agente de cambio y amenaza la diversidad biológica nativa (UICN, 2000).

1.1 ¿Por qué son invasoras?

Existen diversas hipótesis que tratan de explicar los motivos por los cuales algunas especies se convierten en invasoras. La hipótesis mayoritaria es que la falta de enemigos naturales en el nuevo hábitat dónde la invasora llega, le permite crecer sin control invirtiendo menos en compuestos y estrategias defensivas. Esta es la Enemy Release Hypothesis (Catford, Jansson, and Nilsson (2009)).

Otros autores postulan la regla del 10. De diez especies exóticas que llegan a un nuevo terreno, una se convierte en invasora. Esta hipótesis se basa en el azar y el hecho de que algunas comunidades son mas fáciles de ser invadidas que otras (Richardson and Pyšek (2006)).

Una tercera corriente postula que hay ciertos rasgos que proporcionan a las especies invasoras mayores posibilidades para serlo. Por ejemplo, grandes inversiones en dispersión de semillas, ciclos de vida cortos o crecimiento clonal son algunos de los rasgos propuestos para definir un “buen invasor” (Van Kleunen et al. (2010)).

Ninguna de estas hipótesis es excluyente, por lo que los factores determinantes del éxito invasor podría ser una combinación de aspectos de las tres corrientes (Catford, Jansson, and Nilsson (2009)).

Pese a la divergencia de hipótesis que explican la expansión de las especies invasoras, encontramos una mayor unificación con respecto a las teorías que explican el establecimiento de las especies fuera de su rango nativo. No todas las especies pueden ser encontradas en todos los hábitats. Cada especie tiene sus requerimientos y por lo tanto es posible definir su nicho ecológico en base a las condiciones donde una especie puede vivir, crecer y reproducirse (Grinnell (1917)). El nicho de una especie se define como una región (un hipervolumen de n dimensiones) en un espacio multidimensional basado en factores ambientales que afectan al bienestar de la especie (Hutchinson 1957). La similitud climática entre diferentes países podría ser clave para entender el establecimiento de las especies fuera de su rango nativo y una herramienta muy útil para la valoración del riesgo invasor. Un mejor conocimiento de los factores que determinan si una especie podría convertirse en invasora, sería vital para las políticas de conservación de la biodiversidad, pues la prevención es la acción más económica para erradicar una especie invasora (Kumschick and Richardson (2013)). Saber hasta qué punto la similitud con el marco climático de una región, en nuestro caso la Península Íbérica, puede favorecer a la invasión biológica nos puede permitir definir medidas de prevención mayores.

2 Objetivos del trabajo y definiciones iniciales

2.1 Objetivos

El objetivo principal de este trabajo es utilizar diferentes herramientas de gestión y análisis de datos para valorar el riesgo invasor en la Península Ibérica, en base a la hipótesis de la similitud climática como principal factor determinante para el establecimiento de espécies exóticas. Los objetivos específicos son:

Definir el marco climático delimitado por el espacio geográfico de la Península Ibérica.
Cuantificar el grado de similitud climática de todos los países del mundo con el clima de la Península Ibérica.
Determinar si los países con mayor similitud climática son tambien fuente de mayores espécies invasoras en España.

2.2 Definiciones

Marco climatico Ibérico: Espacio multidimensional definido por los valores de diferentes variables bioclimáticas georreferenciados dentro de la Península Ibérica.

Variable bioclimática: Variable descriptora ambiental considerada determinante para el crecimiento, desarrollo y reproducción de la biosfera.

Nicho climático: Espacio multidimensional definido por los valores de diferentes variables bioclimáticas de los puntos de presencia de una especie.

3 Paquetes necesarios

La versión de R con la que se ha trabajado es:

Rv <- R.Version()
Rv$version.string

## [1] "R version 3.5.0 (2018-04-23)"

Antes de empezar, instalaremos los paquetes necesarios a lo largo del documento.

list.of.packages <- c("DataExplorer", "raster", "dismo", "plotly", "fmsb", "factoextra", "MASS","hypervolume","alphahull", "dplyr", "rgdal", "leaflet", "RColorBrewer", "devtools")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

Una vez disponibles todos los paquetes, podemos cargar las respectivas librerías.

library(raster)
library(sp)
library(dismo)
library(ggplot2)
library(plotly)
library(MASS)
require(MASS)
library(fmsb)
library(devtools)
install_github("ggbiplot","vqv")
library(ggbiplot)
library(factoextra)
library(readr)
library(hypervolume)
library(alphahull)
library(DataExplorer)

4 Delimitación del MCI

En primer lugar descargaremos el dataset de la base de datos globales WorldClim que ofrece 19 variables bioclimáticas con una resolución geográfica de 10’ (minutos de grado) ya que resoluciones mayores conllevan demasiado esfuerzo computacional. El argumento var = "bio" de la función getData() reclama dichas variables bioclimáticas. Las variables referentes a la temperatura están multiplicadas por 10.

worldata <- getData("worldclim",var="bio", res=10)
worldata_25 <-getData("worldclim",var="bio",res=2.5)

plot (worldata_25$bio1,
   main = "Mapa de la temperatura media (BIO1) global",
   xlab = "Longitud (ºC)",
   ylab = "Latitud (ºC)")

Los datos de worldata se encuentran en formato rasterizado. Ahora debemos seleccionar la extensión correspondiente a la Península Ibérica. Pasaremos la ventana de latitudes y longitudes que enmarcan la Península Ibérica a la función extent(): (-10º W, 4ºE) y (36ºN, 44ºN). Seguidamente, recortaremos la extensión e sobre los datos worldata utilizando la función crop().

e <- extent(-10, 4, 36, 44)
s.crop <- crop(worldata_25, e)

plot (s.crop$bio1,
   main = "Mapa de la temperatura media (BIO1) en la Península Ibérica",
   xlab = "Longitud (ºC)",
   ylab = "Latitud (ºC)")

El archivo s.crop también se encuentra en formato rasterizado. Para trabajar más cómodamente, construiremos un dataframe llamado marcoclim, definiendo los nombres de las variables bioclimáticas:

s.crop <- crop(worldata, e)
marcoclim <- as.data.frame(s.crop)

BIOs <- paste('BIO', (1:19), sep = '')
names(marcoclim) <- BIOs

Finalmente, exportamos el dataframe marcoclim en formato *.csv.

write.csv(marcoclim, file="Datos/MarcoClim.csv")

5 Análisis exploratorio

Antes de continuar trabajando, necesitamos hacernos una idea de cómo son los datos, para ello realizamos un pequeño análisis exploratorio y data wrangling.

5.1 Data wrangling

En la mayoría de ocasiones los datos con los que se trabaja necesitan de una limpieza y una preparación previa para ser analizados. En nuestro caso, eliminamos los missing values del marcoclim.

marcoclim <- marcoclim[complete.cases(marcoclim),] # Elimina NA's.
marcoclim <- na.omit(marcoclim)
head(marcoclim)

##    BIO1 BIO2 BIO3 BIO4 BIO5 BIO6 BIO7 BIO8 BIO9 BIO10 BIO11 BIO12 BIO13
## 52  134   86   40 4578  247   35  212   85  175   193    76  1287   157
## 53  132   89   40 4698  248   30  218   81  174   192    72  1259   153
## 54  129   94   41 4862  251   23  228   76  173   191    67  1212   147
## 55  129   99   41 5005  255   19  236   73  192   192    64  1145   137
## 56  129  103   42 5164  260   16  244   71  194   194    62  1073   127
## 57  127  107   42 5311  264   11  253   68  195   195    59  1015   118
##    BIO14 BIO15 BIO16 BIO17 BIO18 BIO19
## 52    64    23   422   238   274   394
## 53    63    23   412   235   268   385
## 54    61    22   394   229   258   370
## 55    58    21   368   220   220   346
## 56    55    20   341   209   209   320
## 57    54    19   317   202   202   298

5.2 Estructura de los datos

summary(marcoclim) #resumen de cada variable

##       BIO1            BIO2            BIO3            BIO4     
##  Min.   : 24.0   Min.   : 54.0   Min.   :26.00   Min.   :2680  
##  1st Qu.:117.0   1st Qu.: 91.0   1st Qu.:37.00   1st Qu.:5178  
##  Median :136.0   Median :102.0   Median :38.00   Median :5746  
##  Mean   :134.6   Mean   :100.5   Mean   :37.79   Mean   :5632  
##  3rd Qu.:157.0   3rd Qu.:112.0   3rd Qu.:39.00   3rd Qu.:6236  
##  Max.   :185.0   Max.   :129.0   Max.   :47.00   Max.   :7222  
##       BIO5            BIO6             BIO7          BIO8      
##  Min.   :150.0   Min.   :-76.00   Min.   :127   Min.   : -1.0  
##  1st Qu.:258.0   1st Qu.: -1.00   1st Qu.:237   1st Qu.: 83.0  
##  Median :285.0   Median : 17.00   Median :266   Median :106.0  
##  Mean   :284.2   Mean   : 21.17   Mean   :263   Mean   :105.8  
##  3rd Qu.:312.0   3rd Qu.: 46.00   3rd Qu.:293   3rd Qu.:126.0  
##  Max.   :361.0   Max.   : 87.00   Max.   :334   Max.   :203.0  
##       BIO9           BIO10           BIO11            BIO12       
##  Min.   :-20.0   Min.   : 94.0   Min.   :-38.00   Min.   : 221.0  
##  1st Qu.:181.0   1st Qu.:189.0   1st Qu.: 45.00   1st Qu.: 478.0  
##  Median :205.0   Median :210.0   Median : 64.00   Median : 593.0  
##  Mean   :193.2   Mean   :209.6   Mean   : 66.02   Mean   : 670.0  
##  3rd Qu.:232.0   3rd Qu.:234.0   3rd Qu.: 90.00   3rd Qu.: 800.8  
##  Max.   :275.0   Max.   :277.0   Max.   :129.00   Max.   :1522.0  
##      BIO13            BIO14           BIO15           BIO16      
##  Min.   : 31.00   Min.   : 0.00   Min.   :11.00   Min.   : 85.0  
##  1st Qu.: 62.00   1st Qu.: 5.00   1st Qu.:27.00   1st Qu.:164.0  
##  Median : 80.00   Median :14.00   Median :37.00   Median :217.0  
##  Mean   : 88.25   Mean   :19.62   Mean   :39.13   Mean   :240.3  
##  3rd Qu.:106.00   3rd Qu.:30.00   3rd Qu.:52.00   3rd Qu.:289.0  
##  Max.   :236.00   Max.   :90.00   Max.   :77.00   Max.   :621.0  
##      BIO17            BIO18            BIO19      
##  Min.   :  9.00   Min.   : 16.00   Min.   : 66.0  
##  1st Qu.: 36.00   1st Qu.: 39.00   1st Qu.:125.0  
##  Median : 67.50   Median : 76.00   Median :191.0  
##  Mean   : 82.89   Mean   : 91.43   Mean   :210.4  
##  3rd Qu.:117.00   3rd Qu.:130.00   3rd Qu.:262.8  
##  Max.   :304.00   Max.   :304.00   Max.   :621.0

plot_str(marcoclim) #estructura de los datos
plot_histogram(marcoclim)

5.3 Correlación entre variables

cor(marcoclim) # Idea de la correlación de unas variables con otras.

##               BIO1         BIO2         BIO3         BIO4       BIO5
## BIO1   1.000000000 -0.003218943  0.050613531 -0.058277439  0.7212966
## BIO2  -0.003218943  1.000000000  0.168482451  0.770769526  0.6318544
## BIO3   0.050613531  0.168482451  1.000000000 -0.477637743 -0.1019689
## BIO4  -0.058277439  0.770769526 -0.477637743  1.000000000  0.5984012
## BIO5   0.721296606  0.631854374 -0.101968935  0.598401199  1.0000000
## BIO6   0.883805142 -0.407885261  0.167118931 -0.492863922  0.3452915
## BIO7  -0.017527499  0.918590949 -0.228272638  0.954158045  0.6709731
## BIO8   0.456888738 -0.108221983 -0.127929509  0.001236178  0.2364828
## BIO9   0.652759550  0.273745339  0.130033272  0.116211839  0.6607126
## BIO10  0.919547558  0.297076397 -0.144429394  0.336239959  0.9167325
## BIO11  0.942102524 -0.264616419  0.191978714 -0.384576525  0.4689499
## BIO12 -0.348841653 -0.559844960  0.182554789 -0.609402384 -0.6593703
## BIO13 -0.100078545 -0.603344167  0.198089465 -0.670610900 -0.5041127
## BIO14 -0.758129653 -0.351421682 -0.020090498 -0.262785988 -0.8207783
## BIO15  0.762416555 -0.015640429  0.099351207 -0.121696887  0.5703048
## BIO16 -0.072684836 -0.575272942  0.231543534 -0.668608154 -0.4659488
## BIO17 -0.731581881 -0.375530061 -0.004061334 -0.293794281 -0.8213257
## BIO18 -0.687875511 -0.448629819 -0.051967481 -0.326059433 -0.8322554
## BIO19 -0.001603785 -0.467961092  0.278680322 -0.608561397 -0.3417890
...

plot_correlation(marcoclim) #para visualizar mejor la correlación

5.4 Visualización de la exploración

Ahora, podemos representar los datos obtenidos, dependiendo de las variables que nos interesa mostrar:

Para generar un gráfico con temperatura y precipitación:

plot(marcoclim$BIO1~marcoclim$BIO12, pch=19, col="blue", xlab="Precipitación Anual (mm)", ylab="Temperatura anual media (ºC x10)",main="Temperatura y Precipitación en la Península Ibérica")

Podemos graficar las variables teniendo en cuenta su densidad estimada en base a métodos como Kernel Density Estimation. La estimación de la función de densidad para la temperatura media anual seria:

d <- density(marcoclim$BIO1) # devuelve los datos de densidad 
plot(d)

Para generar un gráfico de densidad con dos variables usamos la función kde2d(kernel density two dimensions) del paquete MASS, usando las variables de temperatura media anual y precimitación media anual:

dens <- kde2d(marcoclim$BIO1, marcoclim$BIO12, n=300, lims = c(0,200,0,1400))
image(dens)

filled.contour(dens,
               color.palette=colorRampPalette(c('gray94','blue','yellow','red','darkred')),
               xlab= "Annual Mean Temperature (ºC)",
               ylab= "Annual Precipitation (mm)", main = "Marco Climático Ibérico")

6 Descartar variables muy correlacionadas (opcional)

Al trabajar con 19 variables, algunas producto lineal de otras como hemos visto en el gráfico de correlación, estamos trabajando con información redundante. Podemos basarnos en el criterio VIF (factores de inflación de varianza), un enfoque simple para identificar la colinealidad entre las variables explicativas para descartar algunas y probar el modelo más adelante con todas las variables y con solo las seleccionadas.

La cuantificación de los factores de inflación de la varianza y el descarte paso a paso se pueden realizar con una función descrita en el blog beckmw:

vif_func<-function(in_frame,thresh=10,trace=T,...){
  
  library(fmsb)
  
  if(any(!'data.frame' %in% class(in_frame))) in_frame<-data.frame(in_frame)
  
  #get initial vif value for all comparisons of variables
  vif_init<-NULL
  var_names <- names(in_frame)
  for(val in var_names){
    regressors <- var_names[-which(var_names == val)]
    form <- paste(regressors, collapse = '+')
    form_in <- formula(paste(val, '~', form))
    vif_init<-rbind(vif_init, c(val, VIF(lm(form_in, data = in_frame, ...))))
  }
  vif_max<-max(as.numeric(vif_init[,2]), na.rm = TRUE)
  
  if(vif_max < thresh){
    if(trace==T){ #print output of each iteration
      prmatrix(vif_init,collab=c('var','vif'),rowlab=rep('',nrow(vif_init)),quote=F)
      cat('\n')
      cat(paste('All variables have VIF < ', thresh,', max VIF ',round(vif_max,2), sep=''),'\n\n')
    }
    return(var_names)
  }
  else{
    
    in_dat<-in_frame
    
    #backwards selection of explanatory variables, stops when all VIF values are below 'thresh'
    while(vif_max >= thresh){
      
      vif_vals<-NULL
      var_names <- names(in_dat)
      
      for(val in var_names){
        regressors <- var_names[-which(var_names == val)]
        form <- paste(regressors, collapse = '+')
        form_in <- formula(paste(val, '~', form))
        vif_add<-VIF(lm(form_in, data = in_dat, ...))
        vif_vals<-rbind(vif_vals,c(val,vif_add))
      }
      max_row<-which(vif_vals[,2] == max(as.numeric(vif_vals[,2]), na.rm = TRUE))[1]
      
      vif_max<-as.numeric(vif_vals[max_row,2])
      
      if(vif_max<thresh) break
      
      if(trace==T){ #print output of each iteration
        prmatrix(vif_vals,collab=c('var','vif'),rowlab=rep('',nrow(vif_vals)),quote=F)
        cat('\n')
        cat('removed: ',vif_vals[max_row,1],vif_max,'\n\n')
        flush.console()
      }
      
      in_dat<-in_dat[,!names(in_dat) %in% vif_vals[max_row,1]]
      
    }
    
    return(names(in_dat))
    
  }
  
}

Aplicamos la función a nuestros datos. Ésta va calculando los VIF, y en cada cálculo quita el que tiene mayor valor hasta quedarse con valores de VIF por debajo del umbral definido (en nuestro caso hasta quedarnos con valores inferiores a 5, un umbral recomendado).

vif_func(in_frame=marcoclim,thresh=5,trace=T)

##  var   vif             
##  BIO1  813.522988309366
##  BIO2  239.44782413085 
##  BIO3  34.7762393382378
##  BIO4  1092.17479361358
##  BIO5  Inf             
##  BIO6  Inf             
##  BIO7  Inf             
##  BIO8  3.96369033534115
##  BIO9  4.21583510546391
##  BIO10 3011.52623004908
##  BIO11 1663.61726148902
##  BIO12 236.137816470082
##  BIO13 99.181388951579 
##  BIO14 105.495934592052
##  BIO15 25.7859393870437
##  BIO16 324.243130786138
##  BIO17 172.888788198011
##  BIO18 119.650658542968
##  BIO19 168.424151605773
## 
## removed:  BIO5 Inf 
...

## [1] "BIO2"  "BIO3"  "BIO8"  "BIO9"  "BIO13" "BIO14"

Estas son las variables con las que nos quedariamos para reducir la colinealidad.

7 Reducción de dimensionalidad mediante PCA

Para describir de forma reducida el MCI realizamos una reducción de dimensiones mediante análisis de componentes principales.

7.1 PCA con todas las variables

Realizamos la PCA con nuestros datos utilizando la función prcomp que R tiene incorporada.

pcaclim <- prcomp(marcoclim, center = TRUE, scale. = TRUE) 
print(pcaclim)

## Standard deviations (1, .., p=19):
##  [1] 2.932734e+00 2.389085e+00 1.445999e+00 1.073491e+00 8.027324e-01
##  [6] 6.757993e-01 4.704042e-01 2.338774e-01 1.637139e-01 1.226069e-01
## [11] 1.002573e-01 7.772701e-02 6.887137e-02 5.791342e-02 4.710744e-02
## [16] 4.006634e-02 3.086699e-02 1.353002e-02 1.181938e-15
## 
## Rotation (n x k) = (19 x 19):
##               PC1          PC2         PC3          PC4         PC5
## BIO1   0.26057151  0.238012832 -0.13443828  0.036857434  0.24466095
## BIO2   0.17820964 -0.248355819  0.33374275 -0.262346960  0.22369365
## BIO3  -0.02885208  0.127133703  0.15301532 -0.849036098  0.09461437
## BIO4   0.16752742 -0.312197451  0.18058120  0.297613338  0.15748431
## BIO5   0.32210187 -0.013221687  0.13112561  0.082911866  0.26497258
## BIO6   0.14467891  0.355591401 -0.18936869 -0.002411807  0.09563225
## BIO7   0.18968588 -0.293422195  0.27336908  0.080156006  0.17451898
## BIO8   0.12949820  0.003907452 -0.55982629 -0.045122752  0.17893624
## BIO9   0.22447097  0.144032253  0.20457806 -0.053542475  0.34923933
## BIO10  0.31312982  0.103748692 -0.05145848  0.156047032  0.27156093
## BIO11  0.18679938  0.324969551 -0.18207586 -0.052433109  0.14884430
## BIO12 -0.27141600  0.202084486  0.19080726  0.116922878  0.21182655
## BIO13 -0.20471851  0.293152290  0.19246920  0.149316359  0.07935544
## BIO14 -0.30956271 -0.109145795 -0.10041250  0.001946151  0.31115272
## BIO15  0.20969829  0.279824140  0.16464927  0.099991960 -0.30461907
## BIO16 -0.19302030  0.303187985  0.22790285  0.129339110  0.05195185
## BIO17 -0.31266003 -0.087197203 -0.08955716  0.008136523  0.36966961
## BIO18 -0.30985703 -0.073151928 -0.16220531  0.048298047  0.33412666
## BIO19 -0.15581732  0.310362018  0.31426410  0.099878118  0.06984306
...

plot(pcaclim, type = "l")

summary(pcaclim)

## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5     PC6
## Standard deviation     2.9327 2.3891 1.4460 1.07349 0.80273 0.67580
## Proportion of Variance 0.4527 0.3004 0.1100 0.06065 0.03391 0.02404
## Cumulative Proportion  0.4527 0.7531 0.8631 0.92379 0.95770 0.98174
##                            PC7     PC8     PC9    PC10    PC11    PC12
## Standard deviation     0.47040 0.23388 0.16371 0.12261 0.10026 0.07773
## Proportion of Variance 0.01165 0.00288 0.00141 0.00079 0.00053 0.00032
## Cumulative Proportion  0.99339 0.99626 0.99767 0.99847 0.99899 0.99931
##                           PC13    PC14    PC15    PC16    PC17    PC18
## Standard deviation     0.06887 0.05791 0.04711 0.04007 0.03087 0.01353
## Proportion of Variance 0.00025 0.00018 0.00012 0.00008 0.00005 0.00001
## Cumulative Proportion  0.99956 0.99974 0.99986 0.99994 0.99999 1.00000
##                             PC19
## Standard deviation     1.182e-15
## Proportion of Variance 0.000e+00
## Cumulative Proportion  1.000e+00

Obtenemos la desviación estándar y la proporción de la varianza de cada una de los componentes principales basándose en las 19 variables climáticas. Con PC1 y PC2 se explica un 75,44% de la variabilidad. Teniendo en cuenta el objetivo del trabajo, consideramos adecuado quedarnos con los dos primeros componentes.

Además, podemos calcular el peso de cada variable bioclimática para cada componente principal, con pca$rotation

pcaclim$rotation[,1] # PC1

##        BIO1        BIO2        BIO3        BIO4        BIO5        BIO6 
##  0.26057151  0.17820964 -0.02885208  0.16752742  0.32210187  0.14467891 
##        BIO7        BIO8        BIO9       BIO10       BIO11       BIO12 
##  0.18968588  0.12949820  0.22447097  0.31312982  0.18679938 -0.27141600 
##       BIO13       BIO14       BIO15       BIO16       BIO17       BIO18 
## -0.20471851 -0.30956271  0.20969829 -0.19302030 -0.31266003 -0.30985703 
##       BIO19 
## -0.15581732

pcaclim$rotation[,2] # PC2

##         BIO1         BIO2         BIO3         BIO4         BIO5 
##  0.238012832 -0.248355819  0.127133703 -0.312197451 -0.013221687 
##         BIO6         BIO7         BIO8         BIO9        BIO10 
##  0.355591401 -0.293422195  0.003907452  0.144032253  0.103748692 
##        BIO11        BIO12        BIO13        BIO14        BIO15 
##  0.324969551  0.202084486  0.293152290 -0.109145795  0.279824140 
##        BIO16        BIO17        BIO18        BIO19 
##  0.303187985 -0.087197203 -0.073151928  0.310362018

Para facilitar el entendimiento de la contribución de cada una de las variables en el análisis de componentes principales, Podemos generar plots con las funciones: -fviz_pca_var -ggbiplot

fviz_pca_var(pcaclim,
             col.var = "contrib", # Color by contributions to the PC
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE     # Avoid text overlapping
)

g <- ggbiplot(pcaclim, obs.scale = 1, var.scale = 1, 
               ellipse = TRUE, alpha=0.02, circle = FALSE)
print(g)

7.2 PCA con la función ‘VIF’

También podríamos realizar una PCA solo con las variables que VIF < 5(opcional)

marcoclimVIF <- marcoclim[,c(2,3,8,9,13,15)]

pcaclim2 <- prcomp(marcoclimVIF, center = TRUE, scale. = TRUE) 
print(pcaclim2)

## Standard deviations (1, .., p=6):
## [1] 1.3290614 1.2863099 1.0979518 0.9198756 0.5886380 0.4252514
## 
## Rotation (n x k) = (6 x 6):
##               PC1        PC2         PC3         PC4         PC5
## BIO2   0.59210269 -0.1444876  0.45372154 -0.10819177  0.28852275
## BIO3   0.08306651  0.3467811  0.47623371  0.77388516 -0.02840042
## BIO8   0.16892244 -0.3077534 -0.61995095  0.60118470 -0.01768585
## BIO9   0.50365198  0.4416262 -0.17743188 -0.15289635 -0.70233760
## BIO13 -0.53763457  0.4944007  0.01922602  0.01303984 -0.02774002
## BIO15  0.26695192  0.5698143 -0.38877671 -0.06652759  0.64929780
##               PC6
## BIO2   0.57247155
## BIO3  -0.21526599
## BIO8   0.36149418
## BIO9   0.05624511
## BIO13  0.68206148
## BIO15 -0.16397395

plot(pcaclim2, type = "l")

summary(pcaclim2)

## Importance of components:
##                           PC1    PC2    PC3    PC4     PC5     PC6
## Standard deviation     1.3291 1.2863 1.0980 0.9199 0.58864 0.42525
## Proportion of Variance 0.2944 0.2758 0.2009 0.1410 0.05775 0.03014
## Cumulative Proportion  0.2944 0.5702 0.7711 0.9121 0.96986 1.00000

Como podemos comprobar, usando el criterio VIF la variabilidad explicada en el PCA con dos componentes resulta del 58.05%, un porcentaje inferior que cuando no eliminamos las variables redundantes, por tanto decidimos no trabajar siguiendo este criterio.

8 Obtención de datos de países y proyección sobre el MCI

Llegados a este punto, pretendemos evaluar si las condiciones climáticas de los diferentes países del mundo se parecen a las de la Península Ibérica. Así se podrá determinar después qué especies podrían encontrar condiciones favorables en nuestro país actuando como invasoras y qué especies españolas podrían serlo para otros países. Con estas evidencias se controlarían mejor la entrada y salida de especies invasoras, mejorando las políticas de comercio y transporte por ejemplo. Por eso proyectamos los datos del marco climático de estos países sobre el nuestro.

Preparamos los datos de la Península Ibérica de los dos primeros componentes:

marcoclim <- marcoclim[complete.cases(marcoclim),] # Elimina NA's

e <- pcaclim$x[,1:2]
edf <- as.data.frame(e)
edf$pais <- "PenIberica" #Añadimos una varaible de País, con todas las etiquetes "PenIberica"
head(edf)

##          PC1       PC2       pais
## 52 -5.632903 2.4546673 PenIberica
## 53 -5.482798 2.0945529 PenIberica
## 54 -5.206638 1.5848376 PenIberica
## 55 -4.514175 1.1294438 PenIberica
## 56 -3.987905 0.6762998 PenIberica
## 57 -3.618972 0.1304936 PenIberica

En este blog, encontramos como se realiza la descarga y recorte de los países. Para llamar a un país concreto se usa la nomenclatura de tres letras definida aquí. Hay una base de datos con los nombres y codigos de los 241 países que incluye la lista. Podemos importar esos datos:

CountryCodes <- read_delim("CountryCodes.csv", 
    ";", escape_double = FALSE, trim_ws = TRUE)

## Parsed with column specification:
## cols(
##   Country = col_character(),
##   A2 = col_character(),
##   A3 = col_character(),
##   `Number,` = col_number()
## )

#A3 son los codigos para descargar el país.

Descargamos los límites políticos de un país, por ejemplo Dinamarca, de la siguiente forma:

denmark <- getData('GADM', country='DNK', level=0)
s<- crop(worldata, denmark)
s2 <- mask(s, denmark)
s2<-as.matrix(s2)
marcoclim.denmark <- as.data.frame(s2)
names(marcoclim.denmark) <- BIOs
marcoclim.denmark <- marcoclim.denmark[complete.cases(marcoclim.denmark),]
h<- crop(worldata_25, denmark)
plot (h$bio1,
   main = "Mapa de la temperatura media (BIO1) en Dinamarca",
   xlab = "Longitud (ºC)",
   ylab = "Latitud (ºC)")

Así, obtenemos un dataframe con las variables bioclimáticas de este país.

Proyectamos estos valores sobre el PCA del marco climático ibérico. Generamos la variable país, con el nombre de la etiqueta del mismo. Además combinamos los dos data frames, con los datos de PC1 y PC2 de España y el país seleccionado.

i <- predict (pcaclim, newdata = marcoclim.denmark) 
i <- i[,1:2]

idf <- as.data.frame(i) 
idf$pais <- "Dinamarca" 

df <- rbind(edf,idf)

Generamos un gráfico para visualizar el grado de solapamiento de ambas nubes de puntos usando la función ggplot

p<-ggplot(df, aes(x=PC1, y=PC2, col=pais)) + geom_point() + ggtitle("Proyección de Dinamarca sobre el PCA del Marco Climático Ibérico")
p + labs(x = "PC1 (46.3%)", y = "PC2 (29.1%)" )

Para entender mejor el solapamiento podemos generar gráficos de densidad:

Para el Marco Climático Ibérico:

dense <- kde2d(edf$PC1, edf$PC2, n=300, lims = c(-12,7,-10,10))
image(dense)

filled.contour(dense,
               color.palette=colorRampPalette(c('gray94','blue','yellow','red','darkred')),
               xlab= "PC1",
               ylab= "PC2", main = "Marco climático Ibérico")

Para Dinamarca:

densi <- kde2d(idf$PC1, idf$PC2, n=300, lims = c(-12,7,-10,10))
image(densi)

filled.contour(densi,
               color.palette=colorRampPalette(c('gray94','blue','yellow','red','darkred')),
               xlab= "PC1",
               ylab= "PC2", main = "Marco climático Ibérico del país proyectado")

Pero, ¿Cómo podemos cuantificar el grado de solapamiento, teniendo en cuenta la densidad de los datos?

9 Cuantificación del grado de similitud climática

9.1 Coincidencia climática versus MCI

Con el fin de determinar el grado de solapamiento climático de un país cualquiera respecto al MCI, creamos hipervolúmenes a partir de los dos primeros componentes principales, que explican un 75.44% de la varianza, y los comparamos, obteniendo así el volumen compartido y otros índices de similitud.

Para ello, utilizaremos el paquete hypervolume (Benjamin, Christine, et al. (2017))(Benjamin, Babich, et al. (2017)), un paquete diseñado para estimar la forma y volumen de bases de datos de más de una dimension y realizar operaciones con ellos: intersección, unión, determinación de los componentes únicos, tests de inclusión y detección de agujeros. Usa una aproximación basada en geometria estocástica ([CRAN] (https://cran.r-project.org/web/packages/hypervolume/hypervolume.pdf)).

Con el objetivo de determinar valores de solapamiento del MCI con todos los países del mundo, crearemos un bucle para contrastar país a país. En primer lugar definiremos el hipervolumen definido por los puntos en el PCA del MCI, a lo que denominaremos marco estático. El marco dinámico estará constituido por aquellos elementos en bucle, dinámicos, de los cuales iremos sacando un conjunto de outputs. En las diferentes pestañas definiremos estos elementos.

9.2 Definiendo el marco climático estático

Una vez obtenidos los datos del MCI, realizado el análisis de componentes principales y creado un data frame con los resultados de éste, definiremos el primer hipervolumen. El hipervolumen de n-dimensiones es una medida que generaliza el concepto de volumen a espacios de dimensión. En nuestro caso utilizaremos un Gaussian KDE(kernel density estimation) en el que todos los puntos contribuyen a la densidad de la probabilidad global ya que testeando vemos que es el mejor algoritmo de inclusion de puntos por nuestra tipologia de distribuciones.

En la función hypervolume_gaussian, por defecto, se le asigna un umbral del 95%, asegurando la máxima probabilidad en la estimación de la densidad. Además generaremos un plot para visualizar el hipervolumen en el que se genera un centroide (en rojo al centro del volumen).

hv1sp = hypervolume_gaussian(data=subset(df,pais=="PenIberica")[1:2730, 1:2])

## 
## Building tree... 
## done.
## Ball query... 
## 
## done.
## 
## Building tree... 
## done.
## Ball query... 
## 
## done.
## Requested probability quantile 0.950000, obtained 0.949316 - setting threshold value 0.000001.
##  For a closer match, you can increase num.thresholds in hypervolume_threshold.

plot(hv1sp,show.contour=TRUE,contour.kde.level=0.01)

9.3 Definiendo el marco climático dinámico

El siguiente paso será generar un hipervolumen para cada uno de los países que testearemos frente al MCI. Por ejemplo, creamos un hipervolúmen con los datos de Dinamarca.

hv2dn = hypervolume_gaussian(data=subset(df,pais=="Dinamarca")[1:238, 1:2])

## 
## Building tree... 
## done.
## Ball query... 
## 
## done.
## 
## Building tree... 
## done.
## Ball query... 
## 
## done.
## 
## Building tree... 
## done.
## Ball query... 
## 
## done.
## 
## Building tree... 
## done.
## Ball query... 
## 
## done.
## Requested probability quantile 0.950000, obtained 0.949491 - setting threshold value 0.000183.
##  For a closer match, you can increase num.thresholds in hypervolume_threshold.

plot(hv2dn,show.contour=TRUE,contour.kde.level=0.01)

9.4 Outputs

Una vez obtenidos los dos hipervolúmenes, utilizaremos la función hypervolume_set y hypervolume_overlap_statistics. La primera genera valores estadísticos computando la intersección, la unión y los componentes únicos(diferencia) de los hipervolúmenes, que podemos observar llamando la función get_volume

La segunda función calcula los métricos de superposición para los dos hipervolúmenes. Imprime directamente los índices de similitud entre los dos volúmenes:

Similitud de Jaccard: Volumen de la intersección de un país respeto al otro, dividido por el volumen de la unión de estos países.

Similitud de Sorensen: Dos veces el volumen de la interesección de un país sobre el otro, dividido por la suma del volumen de ambos países.

Fracción Única 1: Volumen del componente único del primer país dividido por el volumen total de éste.

Fracción única 2: Volumen del componente único del segundo país dividido por el volumen total de éste.

Además calcularemos la distancia entre centroides usando la función hypervolume_distance

setH <-hypervolume_set(hv1sp, hv2dn, check.memory = FALSE)

## Choosing num.points.max=25955 (use a larger value for more accuracy.)
## Using minimum density of 159.367907
## Retaining 15092 points in hv1 and 796 points in hv2.
## Beginning ball queries... 
## 
## Building tree... 
## done.
## Ball query... 
## 
## done.
## 
## Building tree... 
## done.
## Ball query... 
## 
## done.
## Finished ball queries.

volumes<-get_volume(setH)
overl<- hypervolume_overlap_statistics(setH)

dist<- hypervolume_distance(hv1sp, hv2dn, type='centroid')

Con los datos estadísticos obtenidos anteriormente, podemos generar una lista para facilitar la visualización de los diferentes índices y volúmenes.

values <- c(volumes, overl, dist)

En esta figura ilustramos los diferentes índices obtenidos:

Métricos de contraste de similitud

9.5 Bucle para la cuantificación sistemática de los hipervolúmenes

En este punto, implementaremos un sencillo for loop para testear la similitud del marco climático de cada país respecto al de la Península Ibérica de forma rápida y sistemática. El MCI será un valor fijo en este bucle, mientras que los inputs seran los datos climáticos de los diferentes países a testear. A partir de este bucle, generaremos como output un DataFrame que recogerá los volúmenes e índices respecto al MCI de la mayoría de países del mundo. Primero cargaremos una lista con los nombres y los acrónimos de todos los países del mundo:

CountryCodes <- read.csv("CountryCodes.csv", header=FALSE, sep=";", skip = 1)

Decimos “mayoría” de países porque los datos climáticos de algunos de ellos no presentan una distribución geográfica suficientemente amplia como para ser procesados correctamente por la función hypervolume_gaussian del paquete hypervolume. Para evitar que el bucle se interrumpa, debemos eliminar manualmente dichos países:

black_list = list('ASM', 'AIA', 'ATA', 'ABW', 'BMU', 'BVT', 'IOT', 'CYM', 'CXR', 'CCK', 'GIB', 'GRD', 'KIR', 'LIE', 'MAC', 'MDV', 'MLT', 'MHL', 'FSM', 'MCO', 'MSR', 'NRU', 'ANT', 'NIU', 'NFK', 'MNP', 'PCN', 'RUS', 'SHN', 'KNA', 'LCA', 'SPM', 'VCT', 'SMR', 'SCG', 'SYC', 'TKL', 'TUV', 'UMI', 'VAT', 'VGB', 'WLF')
CountryCodesSum <- CountryCodes[ ! CountryCodes$V3 %in% black_list, ]

Como puede verse, la mayoría de los países de la lista black_list son muy pequeños, como por ejemplo SMR (San Marino), GIB (Gibraltar) o VAT (Ciudad del Vaticano). También aparecen muchos estados insulares, como es el caso de CYM (Islas Caimán), BMU (Bermudas) o SYC (Seychelles). ATA (Antártida) aparece en esta lista, no por tener una superficie pequeña (que no es el caso), sino por la escasa distribución geogràfica de sus datos climáticos.

Un caso aparte es el de RUS (Rusia), que ha sido colocada en la black_list debido a la extremadamente amplia distrubución geográfica de sus datos climáticos (debido a su gran extensión). Los ordenadores con los que hemos trabajado no disponían de memoria suficiente como para cargar los datos generados por la función hypervolume_gaussian aplicada a Rusia.

A continuación mostraremos el código del blucle que hemos implementado. Es importante aclarar que el siguiente bloque de código ha sido configurado para no ser ejecutado mediante la opción eval = FALSE:

El bucle recorre todos los países en un lapso de unas 7 horas aproximadamente. Esto trabajando con una workstation de la siguientes caracteríticas:

Memoria RAM: 64 GB.
Disco duro SSD: 256 GB.
Procesador: Processador: Intel Core i7-7700 @ 3.6 GHz de 4 núcleos.

Como ya hemos avanzado, el output de nuestro bucle es un DataFrame llamado Volumes que incluye todos los índices y volúmenes proporcionados por la función hypervolume_gaussian. Finalmete, editaremos brevemente el DataFrame Volumes con la intención de llegar a un DataFrame un poco más leíble y friendly.

CountryCodes$V2 <- NULL
CountryCodes$V4 <- NULL
Volumes <- read.csv("Volumes.csv", header = T)
Volumes$X <- NULL
Vol.Names <- c('Vol.MCI', 'Vol.Country', 'Vol.Intersection', 'Vol.Union', 'Vol.Unique.MCI', 'Vol.Unique.Country', 'Jaccard', 'Sorensen', 'NormVol.MCI', 'NormVol.Country', 'Dist.Centroid', 'Country.Acronym')
names(Volumes) <- Vol.Names
rownames(Volumes) <- Volumes$Country.Acronym
names(CountryCodes) = c('Country.Name', 'V3')
Volume_dataset <- merge(Volumes, CountryCodes, by.x = "Country.Acronym", by.y = "V3")
write.csv(Volume_dataset, file = "Volumes_dataset.csv")

10 Determinación del riesgo invasor

10.1 El top 10 de países de mayor riesgo

Los diferentes índices y valores obtenidos anteriormente nos permiten ordenar los países de forma que podamos visualizar aquellos que presentan un mayor grado de similitud climática.

Cada índice nos revela características distintas de la intersección entre los volúmenes. Vamos a explorar los resultados obtenidos considerando el top 10 de países por cada índice.

Volumes_dataset <- read.csv("Volumes_dataset.csv")

rVolCount <- arrange(Volumes_dataset, desc(Vol.Country)) # Ordenado de mayor volumen a menor
rVolCount <- rVolCount$Country.Acronym # Guardamos la columna con los nombres de los países

rintersection <- arrange(Volumes_dataset, desc(Vol.Intersection)) #Ordenado de mayor interesección a menor
rintersection <- rintersection$Country.Acronym

rvoluniqmci <- arrange(Volumes_dataset, Vol.Unique.MCI) # Ordenado de menor a mayor volumen único en la Península Iberica.
rvoluniqmci <- rvoluniqmci$Country.Acronym

rvoluniqcount <- arrange(Volumes_dataset, Vol.Unique.Country) # Ordenado de menor a mayor volumen único en la Península Ibérica.
rvoluniqcount <- rvoluniqcount$Country.Acronym

rjaccard <- arrange(Volumes_dataset, desc(Jaccard)) # Ordenado de mayor a menor grado de similitud de Jacacard
rjaccard <- rjaccard$Country.Acronym

rsorensen <- arrange(Volumes_dataset, desc(Sorensen)) # Ordenado de mayor a menor grado de similitud de Sorensen
rsorensen <- rsorensen$Country.Acronym

rnormvolmci <- arrange(Volumes_dataset, NormVol.MCI) # Ordenado de menor a mayor fracción única de MCI
rnormvolmci <- rnormvolmci$Country.Acronym

rnormvolcount <- arrange(Volumes_dataset, NormVol.Country) # Ordenado de menor a mayor fracción única del país contrastado
rnormvolcount <- rnormvolcount$Country.Acronym

rdistance <- arrange(Volumes_dataset, Dist.Centroid) # Ordenado de menor a mayor distancia entre centroides
rdistance <- rdistance$Country.Acronym

ranking <- data.frame(rVolCount, rintersection, rvoluniqmci, rvoluniqcount, rjaccard, rsorensen, rnormvolmci, rnormvolcount, rdistance)

names(ranking) <- c('Vol.Country', 'Vol.Intersection', 'Vol.Unique.MCI', 'Vol.Unique.Country', 'Jaccard', 'Sorensen', 'NormVol.MCI', 'NormVol.Country', 'Dist.Centroid')

topten <- ranking[1:10,] #Sacamos los 10 primeros países.
print(topten)

##    Vol.Country Vol.Intersection Vol.Unique.MCI Vol.Unique.Country Jaccard
## 1          IND              CHL            CHL                TCA     ESP
## 2          MMR              ESP            ESP                NLD     ARG
## 3          COL              USA            USA                MYT     GRC
## 4          BTN              ARG            ARG                URY     ZAF
## 5          CHL              BOL            BOL                ALA     ITA
## 6          PNG              AUS            AUS                PLW     LBN
## 7          ECU              MEX            MEX                BHR     PRT
## 8          USA              GRC            GRC                AND     ALB
## 9          CRI              ZAF            ZAF                LUX     TUR
## 10         TWN              LBN            LBN                BRB     AUS
##    Sorensen NormVol.MCI NormVol.Country Dist.Centroid
## 1       ESP         CHL             URY           ESP
## 2       ARG         ESP             NLD           ARG
## 3       GRC         USA             ESP           GRC
## 4       ZAF         ARG             PRT           LSO
## 5       ITA         BOL             LSO           PRT
## 6       LBN         AUS             GRC           ZAF
## 7       PRT         MEX             BEL           ALB
## 8       ALB         GRC             ALB           URY
## 9       TUR         ZAF             LUX           ITA
## 10      AUS         LBN             FRA           LBN

10.2 Mapa de riesgo mundial

La mejor forma de presentar los valores obtenidos de similitud entre países es diseñar un mapa mundial con diferentes colores para aquellas zonas con mayor similitud climática con la Península Ibérica, y por lo tanto mayor riesgo invasor.

Para ello, usaremos el paquete leaflet, que nos permite la creación de mapas interactivos. En primer lugar descargamos los datos de fronteras y polígonos de todos los países del mundo desde el portal Web de Thematicmapping.

Necesitaremos también el paquete rgdal, para obtener csv con los datos descargados, con la función readOGR. Los datos deben estar en un repositorio real del ordenador para cargarlos, por ejemplo así:

#world_spdf=readOGR(dsn= "getwd() ",layer="TM_WORLD_BORDERS_SIMPL-0.3")
world_spdf=readOGR(dsn= "C:/WorldBorders",layer="TM_WORLD_BORDERS_SIMPL-0.3")

Combinamos los datos de los polígonos con la variable que queremos visualizar:

ÍNDICE DE JACCARD:

world_spdf@data$number<-c(1:246)

dfJac <- Volumes[,c(7,12)]
names(dfJac) <- c("Jaccard", "ISO3")

world_spdf@data<-merge(world_spdf@data, dfJac, by= "ISO3", all.x= TRUE)
world_spdf@data<-arrange(world_spdf@data, number)
world_spdf@data$Jaccard[ which(world_spdf@data$Jaccard == 0)] = NA

pal <- colorQuantile("YlOrRd", NULL, n = 9)

state_popup <- paste0("<strong> País: </strong>", 
                      world_spdf$NAME, 
                      "<br><strong>Indice de Jaccard: </strong>", 
                      world_spdf$Jaccard)

leaflet(data = world_spdf) %>%
  addProviderTiles("CartoDB.Positron") %>%
  addPolygons(fillColor = ~pal(Jaccard), 
              fillOpacity = 0.8, 
              color = "#BDBDC3", 
              weight = 1, 
              popup = state_popup) %>%
  addLegend(pal= pal, values= ~Jaccard, opacity=0.9, title = "Jaccard (%)", position = "bottomleft" )

El resultado de éste codigo genera un gráfico interactivo que permite visualizar con zoom los valores del índice en los diferentes países. Aquí mostramos una imagen PNG del mismo, pero para acceder al gráfico interactivo seguir este enlace. Los valores NA corresponden a una similitud nula.

Grado de similitud climatica respecto la Península Ibérica

Similtud climática en base al índice de Jaccard

10.3 Origen de las actuales especies invasoras en España

La Ley 42/2007, de 13 de diciembre, del Patrimonio Natural y de la Biodiversidad creó, en su artículo 61.1, el Catálogo Español de Especies Exóticas Invasoras en el que se incluyen todas aquellas especies exóticas invasoras que constituyan, de hecho, o puedan llegar a constituir una amenaza grave para las especies autóctonas, los hábitats o los ecosistemas, la agronomía, o para los recursos económicos asociados al uso del patrimonio natural (Ministerio de Agricultura , Pesca y Alimentación).

El Catálogo Español de Especies Exóticas Invasoras contempla 181 especies de diferentes grupos taxonómicos provinientes de diferentes países.

Con el objetivo de comprobar si aquellos países con mayor coincidencia climática son fuente de más especies actualmente invasoras, se ha realizado una búsqueda en los portales de ISSG y GBIF para definir los países de origen de las actuales invasoras.

El resultado de esta búsqueda es el siguiente data frame:

origeninvasoras <- read_excel("origeninvasoras.xlsx")

Recontamos el numero de ocurrencias de cada país y lo mezclamos con los datos de similitud:

countriesorigen <-origeninvasoras[,7:56]
b<-as.matrix(countriesorigen)
a<-as.vector(b)
a.t<-table(a)
a.df <- as.data.frame(a.t)
df <- mutate_all(a.df, funs(toupper))
names(df) <- c("Country.Name", "NumSpecies")
withsp<-merge(df, Volume_dataset, by="Country.Name", all.y = TRUE)
withsp$NumSpecies <- as.integer(withsp$NumSpecies)

Seleccionamos el top 10 de países que han aportado mayor número de especies invasoras:

ttorigin <- arrange(withsp, desc(NumSpecies))
ttorigin <- data.frame(ttorigin$Country.Acronym, ttorigin$NumSpecies)
ttorigin[1:10,] #Sacamos los 10 primeros países.

##    ttorigin.Country.Acronym ttorigin.NumSpecies
## 1                       USA                  45
## 2                       CHN                  29
## 3                       MEX                  26
## 4                       BRA                  24
## 5                       ARG                  22
## 6                       JPN                  21
## 7                       CAN                  19
## 8                       PRY                  19
## 9                       URY                  18
## 10                      IND                  17

Podemos representar en forma de correlación los resultados obtenidos, para determinar si la mayor similitud climática explica el número de especies invasoras que aporta ese país:

Distáncia entre centroides vs. Número de especies invasoras aportadas

withsp <- withsp[complete.cases(withsp),] 

p <- ggplot(withsp, aes(x=NumSpecies, y=Dist.Centroid)) +
  geom_point() +
  geom_smooth(method =lm, se=FALSE )
p + labs(x = "Número de especies invasoras en España del país considerado", y = "Distáncia entre centroides climáticos" )

Valor de Índice de Jaccard vs. Número de especies invasoras aportadas

p <-ggplot(withsp, aes(x=NumSpecies, y=Jaccard)) +
  geom_point() +
  geom_smooth(method =lm, se=FALSE)
p + labs(x = "Número de especies invasoras en España del país considerado", y = "Índice de Jaccard" )

10.4 Refinando el riesgo invasor en base a otras variables

A parte de la similitud del marco climático, hay otros factores que pueden resultar determinantes para la determinación del riesgo invasor. La cantidad de mercancías y personas que se desplazan de un país a otro, podrían ser clave para entender en que medida realmente podrían ocurrir los desplazamientos de especies exóticas. Además, en los resultados presentados en el apartado anterior, hay una variabilidad muy grande, probablemente por trabajar con datos de especies invasoras de diferentes grupos taxonómicos. Todo esto requiere la obtención de nuevos datos y desarrollo de nuevos modelos. De forma exploratoria presentaremos algunos modelos simples explorando datos disponibles como son la biodiversidad de cada país (que implicaría mayor probabilidad de existir una espécie que pueda ser invasora) y la distáncia entre el país y la Península Ibérica (lo que probablemente facilite también la llegada de especies alóctonas).

BIODIVERSIDAD: Existen datos de biodiversidad de especies normalizados por el tamaño del país: el National Biodiversity Index. Los valores estan normalizados entre 1 (Màximo: Indonesia) y 0 (Mínimo: Groenlándia).

Importamos los datos

CountryData <- read_excel("Datos/CountryDataSum.xlsx")

names(CountryData)[3] <- "Country.Acronym"
complete<- merge(Volume_dataset, CountryData, by = "Country.Acronym")
complete <-merge(withsp, complete, by="Country.Acronym")

Generamos el gráfico:

complete$nbi <- as.numeric(complete$NBI)

## Warning: NAs introducidos por coerción

p <- ggplot(complete, aes(x=NumSpecies, y=nbi)) +
  geom_point() +
geom_smooth(method = loess, se= TRUE)
p + labs(x = "Número de especies invasoras en España del país considerado", y = "Índice de Biodiversidad nacional" )

## Warning: Removed 25 rows containing non-finite values (stat_smooth).

## Warning: Removed 25 rows containing missing values (geom_point).

Los 10 países con mayor valor de biodiversidad son:

nbirank <- arrange(complete, desc(nbi))
nbirank <- data.frame(nbirank$Country.Acronym, nbirank$nbi)
nbirank[1:10,] #Sacamos los 10 primeros países.

##    nbirank.Country.Acronym nbirank.nbi
## 1                      IDN       1.000
## 2                      COL       0.935
## 3                      MEX       0.928
## 4                      BRA       0.877
## 5                      ECU       0.873
## 6                      AUS       0.853
## 7                      VEN       0.850
## 8                      PER       0.843
## 9                      CHN       0.839
## 10                     CRI       0.820

DISTÁNCIA GEOGRÁFICA: A través del portal Web Geodatos hemos podido obtener los valores de distancia geográfica de todos los países respecto la Península Ibérica.

complete$distance.km <- as.numeric(complete$distance.km)

p <- ggplot(complete, aes(x=NumSpecies, y=distance.km)) +
  geom_point() +
geom_smooth(method = lm, se= FALSE)
p + labs(x = "Número de especies invasoras en España del país considerado", y = "Distáncia respecto el país de origen (km)" )

11 Conclusiones

Durante el desarrollo de este trabajo, hemos aprendido a utilizar distintas herramientas computacionales para valorar el grado de similitud climática de la Península Ibérica con el resto del mundo, con el objetivo de estimar el riesgo invasor. Distintos paquetes de Ry algunas bases de datos disponibles on-line, referentes a variables bioclimáticas de las distintas regiones del planeta, nos han permitido desarrollar un estudio que consideramos muy útil para determinar riesgo invasor, y que sienta las bases de posibles proyectos futuros de evaluación de este tipo. Estos resultados, ayudarían al entendimiento de la problemática relacionada con la pérdida de hábitat y diversidad y podrían mejorar políticas de gestión que protejan o controlen la entrada y salida de especies con potencial invasor, protegiendo así la biodiversidad de nuestro país y también del resto.

Fijándonos en los rankings que hemos obtenido, podemos encontrar diferentes países con más riesgo respecto a la Península Ibérica, dependiendo del métrico o índice al que atendamos. En nuestro caso, hemos determinado el índice de Jaccard como el parámetro más representativo de la similitud de volúmenes para la selección de estos países con condiciones más similares y por tanto de riesgo. Como podemos apreciar en el último mapa, lo países ordenados de mayor a menor riesgo para la Península Ibérica son:

Argentina
Grecia
Sudáfrica
Italia
Líbano
Portugal
Albania
Turquía
Australia.

El resultado de este ranking tiene sentido pues todos los países listados, contienen clima mediterráneo, o bien comparten cercanía con España, o latitud geográfica.

Teniendo en cuenta los otros rankings de número de especies fuente, de biodiversidad y de distancia geográfica, podemos ir más allá en la evaluación del riesgo invasor. En primer lugar, la mitad de países que aparecen en la lista de los 10 países que más especies invasoras han aportado a España, también aparecen en los top 10 de los parámetros de similitud medidos. Por ejemplo, Argentina es el quinto país del que provienen más especies invasoras y a la vez es el país con mayor índice de Jaccard, es decir, más parecido climáticamente a la península ibérica. En segundo lugar, atendiendo la biodiversidad presente respecto a la superficie del país, podemos observar como uno de los países que se encuentra en el top 10 del índice de Jaccard, Australia, también lo está en el top 10 de mayor biodiversidad. Es decir, Australia presenta un riesgo añadido a parte de la similitud climática ya que aumenta la probabilidad de que alguna de las especies llegue a ser invasora en la península ibérica. Otro aspecto interesante a comentar en referencia a los datos de biodiversidad, sería que de la lista del top 10 de países con mayor volumen climático (proyectado sobre la península ibérica), algunos se corresponden con los países de mayor biodiversidad mundial, confirmando la relación ya descrita de mayor biodiversidad en lugares climáticamente más diversos.

Si asumimos que el número de especies invasoras en España viene definido por distintos factores, podríamos evaluar la correlación de cada uno de ellos con ese número. Por ejemplo, según nuestros resultados observamos que a mayor similitud (índice de Jaccard y distancia entre centroides), mayor número de especies invasoras provenientes de ese país. A mayor biodiversidad nacional se ha observado una tendencia a incrementar el número de especies invasoras que han llegado a la península ibérica. Por último, atendiendo a la distancia entre otros países con España observamos que a menor distancia, menor número de especies invasoras seguramente porque los países o regiones cercanas forman parte de la distribución ecológica de las especies que se encuentran en nuestro país o dicho de otra forma, no tendrían carácter invasor por su probable cercanía filogenética y similitud biogeográfica.

Finalmente, nos gustaría resaltar las múltiples posibilidades en cuanto a las variables que se podrían tener en cuenta para seguir evaluando el riesgo, teniendo en cuenta lo desarrollado y obtenido en este trabajo. Por ejemplo, testeo del riesgo invasor de especies animales o plantas a nivel individual y taxonómico o la determinación del riesgo en base al movimiento de ciudadanos o mercancía entre países.

Visto en perspectiva, el presente trabajo representa una primera aproximación de un método para valorar el riesgo invasor desde un punto de vista del data science, una ciencia emergente con gran cantidad de herramientas y en constante evolución que con el tiempo, permitirán el desarrollo de modelos más transversales y de mayor complejidad.

Bibliografía

Benjamin, Blonder, Morrow Cecina Babich, Maitner Brian, Harris David J., Lamanna Christine, Violle Cyrille, Enquist Brian J., and Kerkhoff Andrew J. 2017. “New Approaches for Delineating N-Dimensional Hypervolumes.” Methods in Ecology and Evolution 9 (2): 305–19. doi:10.1111/2041-210X.12865.

Benjamin, Blonder, Lamanna Christine, Violle Cyrille, and Enquist Brian J. 2017. “Using N-Dimensional Hypervolumes for Species Distribution Modelling: A Response to Qiao et Al. ().” Global Ecology and Biogeography 26 (9): 1071–5. doi:10.1111/geb.12611.

Cardinale, B.J., J.E. Duffy, A. Gonzalez, D.U. Hooper, C. Perrings, P. Venail, A. Narwani, et al. 2012. “Biodiversity Loss and Its Impact on Humanity.” Nature 486 (7401): 59–67. doi:10.1038/nature11148.

Catford, J.A., R. Jansson, and C. Nilsson. 2009. “Reducing Redundancy in Invasion Ecology by Integrating Hypotheses into a Single Theoretical Framework.” Diversity and Distributions 15 (1): 22–40. doi:10.1111/j.1472-4642.2008.00521.x.

Grinnell, Joseph. 1917. “The Niche-Relationships of the California Thrasher.” The Auk 34 (4). American Ornithological Society: 427–33. http://www.jstor.org/stable/4072271.

Kumschick, S., and D.M. Richardson. 2013. “Species-Based Risk Assessments for Biological Invasions: Advances and Challenges.” Diversity and Distributions 19 (9): 1095–1105. doi:10.1111/ddi.12110.

Richardson, D.M., and P. Pyšek. 2006. “Plant Invasions: Merging the Concepts of Species Invasiveness and Community Invasibility.” Progress in Physical Geography 30 (3): 409–31. doi:10.1191/0309133306pp490pr.

Simberloff, Daniel, Jean-Louis Martin, Piero Genovesi, Virginie Maris, David A. Wardle, James Aronson, Franck Courchamp, et al. 2013. “Impacts of Biological Invasions: What’s What and the Way Forward.” Trends in Ecology & Evolution 28 (1): 58–66. doi:https://doi.org/10.1016/j.tree.2012.07.013.

Van Kleunen, M., W. Dawson, D. Schlaepfer, J.M. Jeschke, and M. Fischer. 2010. “Are Invaders Different? A Conceptual Framework of Comparative Approaches for Assessing Determinants of Invasiveness.” Ecology Letters 13 (8): 947–58. doi:10.1111/j.1461-0248.2010.01503.x.

Van Wilgen, B.W., S.J. Davies, and D.M. Richardson. 2014. “Invasion Science for Society: A Decade of Contributions from the Centre for Invasion Biology.” South African Journal of Science 110 (7-8). doi:10.1590/sajs.2014/a0074.

Valoración del riesgo invasor en la Península Ibérica en base a la similitud climática

Gianluca Arauz, Olga Burgos, Ramon Diaz y Erola Fenollosa

July, 2018